home *** CD-ROM | disk | FTP | other *** search
/ 8bitfiles.net/archives / archives.tar / archives / commodore-users-of-norman / CUON_121_(06-1985).d64 / fourier smooth (.txt) < prev    next >
Commodore BASIC  |  2019-04-13  |  3KB  |  162 lines
  1. 10 REM - FROM: 'FOURIER SMOOTHING WITHOUT THE FAST FOURIER TRANSFORM',
  2. 20 REM -   BY ERIC E. AUBANEL & KEITH B. OLDHAM ; 'BTYE', FEB 1985, P. 207
  3. 30 REM
  4. 40 REM - MODIFIED FOR THE COMMODORE-64 BY J.J. WORMSER, APRIL 1985
  5. 100 PRINT"[147]"
  6. 110 PRINT"    FOURIER SMOOTHING PROGRAM "
  7. 120 PRINT:PRINT" THIS PROGRAM CALCULATES THE OPTIMUM"
  8. 130 PRINT" SMOOTH CURVE USING THE FOURIER INVERSE"
  9. 140 PRINT" TRANSFORM IN A QUADRATIC DIGITAL"
  10. 150 PRINT" FILTER PROCESS. THE PROGRAM ALSO"
  11. 160 PRINT" INCLUDES A STRAIGHT-LINE CORRECTION"
  12. 170 PRINT" ROUTINE WHICH REDUCES 'END EFFECT'.":PRINT
  13. 180 PRINT" YOU SIMPLY ENTER THE DATA AS PROMPTED. "
  14. 190 PRINT" THE NUMBER OF TRANSFORM POINTS TO BE"
  15. 200 PRINT" KEPT ('E') DETERMINES THE DEGREE OF"
  16. 210 PRINT" SMOOTHNESS. THIS NUMBER MUST LIE"
  17. 220 PRINT" BETWEEN 1 AND INT((N+1)/2). THE"
  18. 222 PRINT" DEGREE OF SMOOTHNESS IS INVERSELY"
  19. 224 PRINT" PROPORTIONAL TO THIS VALUE."
  20. 226 PRINT
  21. 230 PRINT" HIT ANY KEY TO CONTINUE";
  22. 240 GETK$:IFK$=""THEN240
  23. 250 PRINT"[147]":CLR
  24. 260 INPUT" ENTER NUMBER OF DATA POINTS";N
  25. 270 REM LEAVING R&I ARRAYS UNDIMENSIONED    LIMITS VALID VALUES OF E TO <=10
  26. 280 N2=INT((N+1)/2+1):DIMX(N),X1(N),U(N2),V(N2)
  27. 290 FOR I=0TON-1
  28. 300 INPUT" ENTER DATAPOINT VALUE";X(I)
  29. 310 NEXTI
  30. 320 GOSUB450
  31. 330 PRINT"[147]":PRINT" SMOOTHED DATA FOR 'E'=";ES;"ARE:":PRINT
  32. 340 FORI=0TON-1
  33. 350 PRINT" X(";I;")=";X1(I)
  34. 360 NEXTI
  35. 370 PRINT:PRINT" IF YOU WANT TO TRY A DIFFERENT E,"
  36. 380 INPUT" ENTER (1) ELSE ENTER (0)";MO
  37. 390 IF MO<1ORMO>1THEN420
  38. 400 IF  MO=1THENGOSUB450
  39. 410 GOTO330
  40. 420 END
  41. 430 REM FOURIER ALGORITHM SUBROUTINE        BEGINS AT LINE 60. LINE #'S ARE THE
  42. 440 REM SAME AS FOR THE HP VERSION OF THE   SUBROUTINE
  43. 450 PI=3.141593
  44. 460 PRINT"[147]":PRINT" NUMBER OF TRANSFORM POINTS 'E'";
  45. 470 INPUT ES
  46. 475 PRINT
  47. 480 IFES>INT((N+1)/2)THENPRINT" 'E' TOO LARGE":FORI=1TO1000:NEXT:GOTO460
  48. 490 IFES<>INT(ES)ORES<=1THENGOTO460
  49. 500 IFES<=QTHEN1260
  50. 510 REM
  51. 520 IFQ=<>0THEN720
  52. 530 REM STRAIGHT-LINE CALCULATION
  53. 540 S1=0:S2=0
  54. 550 D=INT(N/10)
  55. 560 FORJ=0TOD-1
  56. 570 S1=S1+X(J)
  57. 580 S2=S2+X(N-J-1)
  58. 590 NEXTJ
  59. 600 X1=S1/D:X2=S2/D
  60. 610 M=(X2-X1)/(N-D)
  61. 620 B=(X1+X2)/2-M*N/2
  62. 630 REM CALCULATE(R0)
  63. 640 G=0
  64. 650 FORJ=0TON-1
  65. 660 X(J)=X(J)-M*J-B
  66. 670 G=G+X(J)
  67. 680 NEXTJ
  68. 690 R(0)=G/N
  69. 700 Q=1
  70. 710 REM
  71. 720 PRINT" WORKING ON CALCULATIONS. PLEASE WAIT..."
  72. 730 PRINT
  73. 740 J2=INT((N-1)/2)
  74. 750 P1=INT(100*LOG(2*J2-1)/LOG(2)+.5)/100
  75. 760 FORK=QTOES-1
  76. 770 J1=J2
  77. 780 S=PI*K*2/N
  78. 790 C=COS(S):S=SIN(S)
  79. 800 FORJ=1TOJ1
  80. 810 L=2*J-1
  81. 820 U(J)=X(L)*C+X(L+1)
  82. 830 V(J)=X(L)*S
  83. 840 NEXTJ
  84. 850 S=2*S*C:C=2*C*C-1
  85. 860 FORP=1TOP1
  86. 870 U(J1+1)=0:V(J1+1)=0
  87. 880 J1=INT((J1+1)/2)
  88. 890 FORJ=1TOJ1
  89. 900 L=2*J-1
  90. 910 U=U(L)*C-V(L)*S+U(L+1)
  91. 920 V(J)=U(L)*S+V(L)*C+V(L+1)
  92. 930 U(J)=U
  93. 940 NEXTJ
  94. 950 S=2*S*C:C=2*C*C-1
  95. 960 NEXTP
  96. 970 R(K)=(X(0)+(U(1)*C+V(1)*S))/N
  97. 980 NEXTK
  98. 990 REM
  99. 1000 FORK=QTOES-1
  100. 1010 J1=J2
  101. 1020 S=2*PI*K/N
  102. 1030 C=COS(S):S=SIN(S)
  103. 1040 FORJ=1TO J1
  104. 1050 L=2*J-1
  105. 1060 U(J)=-(X(L)*S)
  106. 1070 V(J)=X(L)*C+X(L+1)
  107. 1080 NEXTJ
  108. 1090 S=2*S*C:C=2*C*C-1
  109. 1100 FORP=1TOP1
  110. 1110 U(J1+1)=0:V(J1+1)=0
  111. 1120 J1=INT((J1+1)/2)
  112. 1130 FORJ=1TOJ1
  113. 1140 L=2*J-1
  114. 1150 U=U(L)*C-V(L)*S+U(L+1)
  115. 1160 V(J)=U(L)*S+V(L)*C+V(L+1)
  116. 1170 U(J)=U
  117. 1180 NEXTJ
  118. 1190 S=2*S*C:C=2*C*C-1
  119. 1200 NEXTP
  120. 1210 I(K)=-((U(1)*C+V(1)*S)/N)
  121. 1220 NEXTK
  122. 1230 REM
  123. 1240 IFES>QTHENQ=ES
  124. 1250 REM
  125. 1260 REM
  126. 1270 REM
  127. 1280 REM CALCULATE X1(0)
  128. 1290 F1=0:F2=0
  129. 1300 FORK=1TOES-1
  130. 1310 T=R(K)
  131. 1320 F1=F1+T
  132. 1330 F2=F2+K*K*T
  133. 1340 NEXTK
  134. 1350 X1(0)=R(0)+2*(F1-F2*(1/ES/ES))
  135. 1360 X1(0)=X1(0)+B
  136. 1370 REM
  137. 1380 P1=INT(100*LOG(2*ES-3)/LOG(2)+.5)/100
  138. 1390 FORJ=1TON-1
  139. 1400 T2=ES*ES
  140. 1410 FORK=1TOES-1
  141. 1420 F=1-K*K/T2
  142. 1430 U(K)=R(K)*F:V(K)=-(I(K)*F)
  143. 1440 NEXTK
  144. 1450 K1=ES-1
  145. 1460 S=2*PI*J/N
  146. 1470 C=COS(S):S=SIN(S)
  147. 1480 FORP=1TOP1
  148. 1490 U(K1+1)=0:V(K1+1)=0
  149. 1500 K1=INT((K1+1)/2)
  150. 1510 FORK=1TOK1
  151. 1520 L=2*K-1
  152. 1530 U=U(L)*C-V(L)*S+U(L+1)
  153. 1540 V(K)=U(L)*S+V(L)*C+V(L+1)
  154. 1550 U(K)=U
  155. 1560 NEXTK
  156. 1570 S=2*S*C:C=2*C*C-1
  157. 1580 NEXTP
  158. 1590 X1(J)=R(0)+2*(U(1)*C+V(1)*S)
  159. 1600 X1(J)=X1(J)+M*J+B
  160. 1610 NEXTJ
  161. 1620 RETURN
  162.